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ABSTRACT 

Dust grains migrating under Poynting-Robertson drag may be trapped in mean-motion 
resonances with planets. Such resonantly trapped grains are observed in the solar sys¬ 
tem. In extrasolar systems, the exozodiacal light produced by dust grains is expected 
to be a major obstacle to future missions attempting to directly image terrestrial plan¬ 
ets. The patterns made by resonantly trapped dust, however, can be used to infer the 
presence of planets, and the properties of those planets, if the capture and evolution of 
the grains can be modelled. This has been done with N-body methods, but such meth¬ 
ods are computationally expensive, limiting their usefulness when considering large, 
slowly evolving grains, and for extrasolar systems with unknown planets and parent 
bodies, where the possible parameter space for investigation is large. In this work, 
we present a semi-analytic method for calculating the capture and evolution of dust 
grains in resonance, which can be orders of magnitude faster than N-body methods. 
We calibrate the model against N-body simulations, finding excellent agreement for 
Earth to Neptune mass planets, for a variety of grain sizes, initial eccentricities, and 
initial semimajor axes. We then apply the model to observations of dust resonantly 
trapped by the Earth. We find that resonantly trapped, asteroidally produced grains 
naturally produce the ‘trailing blob’ structure in the zodiacal cloud, while to match 
the intensity of the blob, most of the cloud must be composed of cometary grains, 
which owing to their high eccentricity are not captured, but produce a smooth disk. 


1 INTRODUCTION 

Small dust grains are released into the Solar system in colli¬ 
sions between asteroids and by the outgassing of comets. If 
they are not so small as to be unbound, they will spiral to¬ 
wards the Sun due to Poynting-Robertson (PR) drag. These 
dust grains reflect, absorb, and re-radiate the solar light, 
producing the zodiacal light and zodiacal thermal emission. 
Around other stars, similar processes may produce exozodi¬ 
acal light. Such exozodiacal light is expected to be a signif¬ 
icant barrier to future missions attempting to directly im¬ 
age terrestrial planets if the level of dust is greater than 



aboard New Horizons may soon detect dust in resonance 


with Neptune (Vitense, Krivov & Lohne 2014). As similar 
processes might be expected to occur in extrasolar planetary 
systems, this raises the interesting possibility that structures 
in extrasolar debris disks are created by this phenomenon, 
and can be used to infer the presence of planets and their 
properties. 

In the outer part (>> 10 au) of a stellar system, de¬ 
bris disks can be spatially resolved by facilities and instru- 
ment s such as HST (|Kalas, Graham &; Clampin|2005 ), Sub¬ 


aru (Thalmann et al. 2011), SMA (Hughes et al. 2011), 


VLT (Buenzli et al. 2010), Herschel (Duchene et al. 


Matthews et al. 2014), Gemini (Wahhaj et al. 2014), and 
ALMA ( Boley et al.|2012 ). This could allow the inference of 


2014 


resonantly trapped dust from the spatial distribution (Kuch- 
ner &; Holman||2003 Wyatt [2003 ), although in these bright 
disks collisions are far more important than PR drag (Wyatt 


2005), which may wash out such patterns (Kuchner &; Stark 


2010). Possible detections of clumpy structures which may 
be associated with planetary resonances have already been 
reported ( Holland et al.|[l998} Greaves et al. 2005 Corder 
et al.|2009 ). In the inner parts (<<10 au) of systems, bright 
dust disks are much rarer ( Fujiwara et al.||2013 ), and tran¬ 
sient in nature (Wyatt et al. 2007). However, bright outer 



































































































2 Shannon et al. 


disks may leak dust inwards due to PR drag (Mennesson 
et al.|in p ress), where resonant capture may occur. PR drag 


can also be relatively more important in dimmer disks. The 
Keck Nulling Interferometer (Millan-Gabet et al. pol l 


Men¬ 


nesson et al. 2013| and the Large Binocular Telescope In¬ 
terferometer Hinz 2013; K ennedy et al.|2 014) can allow for 
the detection of these dimmer disks. These interferometers 
detect the flux leak from a dust disk through an interference 
pattern. With repeated measurements of the same system, 
rotating structures, such as resonantly trapped dust, may 
be able to be inferred. 

To model resonantly trapped dust, various groups have 
used N-body simulations of dust grains. |Dermott et al.| 
(1994) simulated the inspiral of 912 dust grains with diam¬ 
eters of 12 pm from the asteroid belt, and showed that the 
pattern produced by the capture and persistence of those 
grains in resonance is compatible with that observed by the 
Infrared Astronomical Satellite. Deller &; Maddison (2005) 
integrated 500 dust grains under PR drag for 120 differ¬ 
ent systems, generating a catalogue of images for different 
planet mass, planet eccentricity, dust size, dust parent body 
eccentricity, and dust parent body semimajor axis. Stark 


&; Kuchner (2008) required a 420 core cluster to perform 


simulations of 5000 particles spiralling past a planet under 
PR drag for 120 different parameter combinations. Similarly, 


Rodigas, Malhotra & Hinz (20141) used a suite of 160 sim¬ 


ulations to produce a fitting formula for disk appearances. 
In these cases, the long computation times of N-body sim¬ 
ulations makes their use to generate potential disk images 
on a system-by-system case somewhat impractical, hence 
the aforementioned groups have generated catalogues of disk 
images and fitting formulae to be matched with disks once 
discovered. 

These high computation times are set by the need to 
resolve individual orbital periods of the planet (s) and dust 
grains. However, the evolution under PR drag and the evo¬ 
lution in resonance occur on much longer time-scales than 
the orbital period. Thus, in this paper we develop a model of 
resonant capture and evolution, with which dust grain tra¬ 
jectories can be calculated on the PR drag time-scale, rather 
than the orbital time-scale, resulting in orders of magnitude 
decreases in the computation time. 

In Sj2] we outline the physics of a single dust grain un¬ 
dergoing PR drag around a star where it may be captured 
into mean motion resonance with a planet. In ^3] we perform 
a suite of iV-body simulations, and use them to validate 
and calibrate our model. In © we describe how to use the 
equations of evolution to produce disk images, and compare 
them to the N-body simulations images. In ^5] we compare 
the output of the model to measurements of the structure 
of Earth’s resonantly trapped ring. In ^ 6 ] we discuss other 
possible applications of the model. 


2 SINGLE PARTICLE PHYSICS 

Consider a star of mass ra*, orbited by a dust grain, the or¬ 
bit of which is described by six parameters: semimajor axis 
a, eccentricity e, inclination z, longitude of pericentre (tu), 
ascending node (Q), and mean anomaly (M). The orbital 
evolution of the dust grain will be dictated by the param¬ 
eter /3 pr, the ratio of the force of radiation pressure to the 


gravitational force. If the star is also orbited by a planet of 
mass rrip, and orbital parameters a p , e p , i p , Q p , M p , res¬ 
onant interactions between the planet and dust grain may 
also be important to its orbital evolution. In this section we 
consider these effects analytically. These expectations will 
then be compared to simulations in §3.3| 


2.1 Poynting-Robertson Drag 

Assuming dust grains to be spherical, they experience a force 
from the stellar radiation of the form 


Eradiation — /3pR Eg ra vity 


1 - 2 - 


f - ( ^ 1 Q 


(1) 


where /3pr is the ratio of the force of radiation pressure to 
gravity acting on the grain, which will depend on its size, 
shape and composition. For example, a spherical blackbody 
grain of size s and density p has 

3L* 


PPR = 


( 2 ) 


87 vcpGm*s ’ 

where L* is the stellar luminosity. 

Averaged over an orbit, this causes dust particles to 


evolve in semimajor axis and eccentricity as (Wyatt & Whip- 
ple||1950 ) 

Gm, /3 P r (2 + 3e 2 ) 


da 

dt 


-0.624-^-/3 P r 

kyr m© 


(1 -e 2 )2 

m* lau (2 + 3e 2 ) 

(l-e 2 )t : 


(3) 


de 

dt 


5 Gm* /3 pr 
~ 2~a 


— 1.56kyr ^ PR 


c (1 _e 2 )2 

m* /lauV 
m© V a ) 


(1 — e 2 ) 2 


(4) 


Integrating equation [3] from an initial semimajor axis ao and 
eccentricity eo = 0 , a dust particle will reach the star in time 


tpr — 


ale 


4Gm* /3 pr 


= 0.40 (-2^-) —-2—kyrs. 

Vlau/ m* ppr 


(5) 


Particles with higher eo reach the star more quickly, but the 
problem does not lend itself to a compact analytic form. 

During the orbital evolution, there is a constant of in¬ 
tegration which can be obtained from equations [3] and [4] 
( Wyatt &; Whipple) 1950| 

K = a (l — e 2 ) e - M ( 6 ) 

The value of K is constant over many orbits, but small 
variations occur during an orbit. Differentiating equation 
6 ] with respect to time, substituting in a and e from [Burns] 
Lamy & Soter (1979), and taking the expression to first order 


in e, we find 


1 dK 
K df 


8 /3pr Ukepl e 


cos /, 


5 e c 

and thus a libration in the constant K with a magnitude: 

A K ^ 8 ydpR Ukepl. 

K 5 e c 


(7) 


( 8 ) 
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over the course of an orbit. 


2.2 Mean motion resonance 

A dust grain can be in a k -th order mean motion resonance 
with a planet if their mean motions are close to the ratio 
j : j + k, where j and k are positive integers. From the mean 
longitudes A = M + vj we construct a resonant angle 

P — j + p — (j + k) A + km, (9) 


where A and X p are the mean longitude of the dust grain and 
the planet respectively. For our purposes, vj evolves slowly 
compared to A and A p , so the j : j + k resonance occurs at 

2 

aj:j + k ~ (1 - /3pr) 3 a P’ ( 10 ) 


which differs from the usual expression due to the effects 
of radiation pressure. Here the match is only approximate 
as the km term allows for a small offset in a, and there 
is a small libration of a about the nominal resonance lo¬ 
cation corresponding to the libration in c p. While the dust 
grain is near resonance, the perturbations from the planet 
can be approximated by taking the terms in the disturbing 
function that are first order in eccentricity, and using La¬ 
grange’s equations for the evolution of the orbital elements 
to give ( Murray &; Dermott|1999 ) 

= — 2jC r ae k sin (p, (11) 

j:j+k 


dt 


where 

Cr — 


de 

dt 


= kC r e k 1 simp, 


j:j+k 


(«) + /<(«)) 

na z a p 

-VGrrip ( afd ( a ) + fi (op) 

^(2±fc)’(l-^p R )I a 2 


( 12 ) 


(13) 

(14) 


where n — \JGm* (1 — /3pr) a -3 , a — a p /a ~ 
(1 - /3pr) - 1/3 ((j + k) / j )~ 2/3 , fd (a) is the direct term of 
the disturbing function, and fi (a) is the indirect term of 
the disturbing function. 

Neglecting PR drag, and using equation [To| for the dust 
grain’s semimajor axis, taking the second time derivative of 
ip, neglecting the contribution from X p and zu terms, one 
finds ip evolves as 


Cp — 3j 2 C r ne k sin <p, (15) 

.3 

= 3 . J y/ Grn * (1 - /3pR)a~ 1 C r e k sin <p. (16) 

J + k 

As C r < 0, this is the equation for a pendulum with fre¬ 
quency 

c J 2 m 3j 2 \C r \ne k (17) 

.3 

= 3 3 \JGm * (1 - /3 P R )a~ 1 | C r \ e k . (18) 

3 + k 


Thus, p can librate or circulate, with libration in the case 
that the dust grain is in resonance with the planet. This is 
similar to the case without radiation pressure, but \C r \ con¬ 
tains dependence on /3pr, and a , which also depends on 0pr. 


As noted, for a fixed a p there is a small range of possible 
a for the dust grains for which it is still possible to construct 
a resonant angle that can librate. Defining C r as equation 
13\ and otherwise following the derivation from |Murray~fc] 
Dermott (1999), the maximum libration width of a particle 
in a k = 1 resonance is 


Aa max /16 \C r \ \ 2 / 1 \Cr\\ 2 

a j: j + 1 V 3 n e J { + 27j 2 e 3 n J 

and for a k = 2 resonance the width is 

Aa max _ ± /16 \Cr\_e2\ 5 

a j: j +2 V 3 n J 


J^\Cr\ 

9 je n 
(19) 


( 20 ) 


We do not investigate resonances with higher k , as our N- 
body simulations produce very few captures at k > 1. 


2.3 Resonance capture 


As particles migrate past a planet’s resonances due to PR 
drag, the question of whether they become caught in res¬ 
onance is deterministic. However, as it can depend on the 
relative orbital phases of the planet and dust grain as res¬ 
onance is approached, it can be treated probabilistically 
( Henrard|1982| Quillen|2006 ). Mustill & Wyatt (2011) use a 
simple Hamiltonian model which uses the lowest order term 


of the disturbing function (as we did in [2.2) to calculate 
the capture rates and initial libration widths for massless 
particles that encounter first and second order mean mo¬ 
tion resonances during migration, as well as the eccentricity 
kicks during resonant crossings that do not result in cap¬ 
ture. While the model is developed in the context of the 
capture of planetesimals by a migrating planet, the calcu¬ 
lations depend only on the rate of change of the bodies’ 
orbital separation, and so can be applied to the capture of 
migrating dust particles by a fixed planet. However, that 
model assumes /3 pr = 0 , which requires a slight correction 
for our context. This is because when /3 pr > 0, the reduc¬ 
tion in effective stellar mass seen by the dust grains due to 
radiation pressure means that resonances occur closer to the 
planet than when particles do not feel radiation forces, and 
therefore the strength of a given resonance is different for 
particles of different /3 pr. 

We present details of the required changes to the Hamil¬ 
tonian model, including changes to the resonant strengths 
and eccentricity damping effects, in the Appendix. Based on 
these changes, we now predict capture probabilities, as well 
as initial libration widths and eccentricity kicks during non¬ 
capture resonance crossings, for real scenarios. We adopt a 
stellar mass of 1M®, planet mass of Ira©, and planet semi¬ 
major axis of 1 au. We integrate trajectories of dust particles 
with /3 pr = 0.005,0.01,0.02,0.04,0.08,0.16 and 0.32, with 
initial eccentricities of 0.01. Encounters with resonances 
from the 2:1 up to the 19:18 are integrated, unless said reso¬ 
nance lies inside the planet’s orbit. In figure [l] we show the 
capture probabilities for first-order resonances as a function 
of dimensionless parameters J (oc e 2 ) and dB/dt (oc da/c?t)Q 
In this figure we overplot the parameters corresponding to 
capture in the first-order resonances of a 1M© planet for 


1 We use B for the parameter denoted (3 in MW11, to avoid 
confusion. 
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azimuthal velocity (vo), the change in specific energy of the 
dust grain is 


SE 


2 Gm* 
e 2 a p 





(23) 


For an initial ei and a hnal e2, the specific orbital energy 
change is also 


SE 


Gm * (1 - /?) 
2a p 


(ei - £ 2 ), 


(24) 


and equating the two energies yields 

^(wi-frR(i-0) 3 > 

(25) 


Cl — C2 = 


which assumes e± — 62 << e. 

After one conjunction, the planet and dust grain move 
apart in longitude, meaning that the next conjunction will 
happen at 


A c 


27T ^1 - y/l-fr 



-1 


(26) 


Figure 1. Capture probabilities calculated by integrating the 
Hamiltonian model (white = 100%, black = 0%), including eccen¬ 
tricity damping at k = 5 X 10 —5 . The x-axis measures a dimen¬ 
sionless “eccentricity” (J oc e 2 ), and the y-axis a dimensionless 
migration rate ( dB/dt oc da/dt). Overplotted are the dimension¬ 
less momentum and migration rate corresponding to the compar¬ 
ison integrations described in the text. 


dust particles with /3 pr = 0.005 up to (3pr = 0.32. We can 
see that many of these resonances lie in the regime of pos¬ 
sible but not certain capture, meaning that a population of 
drifting dust grains will populate a range of resonances. 


If we assume that interactions between the dust grain 
and the planet will be dephased if the change in the longi¬ 
tude of conjunction SX C is greater than some amount A, this 
yields a constraint that 


12t r 


v/1 — /^PR 




1 - V 1 “ ( 1 “ 2 e 


^ A, 


(27) 


which reduces to the scaling of Wisdom ( |1980 ) in the limit 
/3 pr —>■ 0 of 


2.4 Maximum j of capture 


Wisdom (1980) showed that first order resonances over¬ 
lap when they are at semimajor axes less than ~ 
1.3 (m p /m*) 2 / 7 a p greater than that of the planet. Bodies 
move chaotically between the overlapping resonances, which 
makes them unstable. Consequently, we would expect such 
resonances to not capture migrating dust grains. Here we 


employ a different approach to that of Wisdom (1980) to 


calculate the width of this chaotic zone, but which recovers 
their result, at the same time accounting for the effects of 
radiation pressure. 

In the inertial frame of a planet with semimajor axis a p , 
a dust particle with radiation pressure co-efficient /3 pr on 
a circular orbit with semimajor axis a p (1 + e) approaches 
from infinity with velocity 


V °° “ ^|^(Wl-0PR ( X - I)) ’ ( 21 ) 

assuming e << 1, an assumption we make throughout this 
derivation. The particle’s path is deflected by an angle 0, 
given by 

9 - ~I)Y 2 < 22 > 

assuming the angle is small. Converting back to the frame 
of the star, and considering the final radial velocity (v r ) and 


e<f^Y f^) f 


[- 3 A J 


(28) 


A reasonable expectation is that A — 2i r, which results in 
a leading coefficient of « 1.5, slightly larger than the value 


obtained by Wisdom (1980), though numerically, Duncan, 


Quinn &; Tremaine (1989) found this coefficient to be ~ 1.5, 
and Chiang et al. (2009) preferred a coefficient of ~ 2.0 to 


fit their simulations of Fomalhaut’s disk. Recent results have 
also shown some dependence on system age and collision 


time (Nesvold & Kuchner 2014 Morrison & Malhotra 2014|. 

For our own simulations, we will fit this coefficient in §3.3.2| 


2.5 Resonant evolution 


In the absence of drag a dust grain resonating with an in¬ 
terior planetary perturber librates around (p = 7r (Equation 
With the addition of PR drag, the centre of libration 
changes. Momentum and energy loss to PR drag must be 
offset by a gain of momentum and energy from the resonant 
interaction. Equating equations [3] and m we can estimate 
that (p librates around (/?o, given by 


po 


sin 


sin 


k ^kepler,j:j-|-k /^PR T 3 e ) \ 

C r c 2je k (1 _ e 2) | J ’ 


(1 — (9pr) 3 /3pr 2 + 3e 2 Gm * j 3 


2 e k C r 


(j + k) 


(29) 
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Note that this expression does not hold for the 2 : 1 
resonance, or other resonances which exhibit asymmetric 
libration for which higher order terms in the disturbing 
function are important (Message 1958). For those reso¬ 


nances there are two possible centres of libration with differ¬ 
ing capture probabilities (Mu rray-Clay &; Chiang||2Q05|. In 
this case, a numerically calibrated prescription from [Wy-] 


att 


(2003) may be employed, which we do in £j4| They 
found that the centre of libration is given by cos^q 1 = 
0.39 — 0.061/e. Using the PR-drag migration rate of equa¬ 
tion [3] their results would predict that the relative chance 
of capture into the lower 2 : 1 resonance is 1/2 — 

0.85yGm*/3 PR /c(m p a) -0 ' 25 (2 + 3e 2 ) 0 ' 5 (l - e 2 ) -0 ' 75 . 

While the dust particle is in resonance, it undergoes 
only small variations in a. The eccentricity, however, may 
evolve substantially. Taking the value of ipo from equation |29| 
and substituting it into equation [12] gives an average e from 
resonant interaction. Combined with the e from PR drag 
(equation [4]) , the eccentricity evolves according to 


de _ 1 Gm* /3pr 
dt 2 a 2 


k (2 + 3e 2 ) 


5e 


% • k c \(j + k)e{ 1 — e 2 )i 


(30) 


Expanding to second order in e and solving yields (Weiden- 
schilling & Jackson 1993) 


e(t) = VWrfc) (31) 

Note that the expression in |Weidenschilling &; Jackson| 
(1993) has a typographical error: it is missing the square 
root. Here the growth time of the particle’s eccentricity is 


= 0.2 a 2 j:j+k c/ (Gm^PR). 


(32) 


motion resonance with the planet. For a particle not cap¬ 
tured in lower j resonances, the method of §2. 3| calculates a 
probability of capture into the 6:5 resonance of ~ 0.3. Allow¬ 
ing for capture into lower j resonances, the overall capture 
probability into the 6:5 mean motion resonance for particles 
with similar initial conditions is ~ 0.24. Thus, we expect this 
is a typical outcome. After capture, the semimajor axis of 
the dust particle oscillates, but remains within the possible 
values for a resonant dust grain as outlined in equation [T9| 
(top left panel of figure [2|. After capture, the evolution of 
the eccentricity is well described by equation [3l] (top right 
panel of figure [2} . The expected centre of libration of the 
resonance angle (c^o) changes quickly at early times, as ex¬ 
pected from equation [29] while the centre of libration of 
the resonance angle in the simulation catches up on a much 
slower timescale (bottom left panel of figure[2]). The libration 
width grows exponentially with time (bottom right panel of 
figure [21. At 3.6 x 10 5 years, the particle escapes from reso¬ 
nance. Later on we will quantify when and how dust grains 
escape from resonance. 

After escaping from resonance, the particle continues 
to evolve under PR drag decreasing its semimajor axis as 
equation [3] and its eccentricity as equation [4] At about 3.9 x 
10 5 years, the particle is removed from the simulation once 
a ~ 0.05, at which point the integration timestep is too large 
to resolve the orbit. 

From this example, we can see that we have some un¬ 
derstanding of the capture process. However, as we have 
noted, some elements are not well understood (libration 
width growth, resonance escape), and others (capture proba¬ 
bility, resonance widths at capture) would benefit from addi¬ 
tional validation. To this end, we perform a suite of N-body 
simulations varying system parameters of interest. 


3 N-BODY INTEGRATIONS 


3.1 Single particle example: 
drag 


Evolution under PR 


To enhance qualitative understanding, we begin by present¬ 
ing the evolution of a single dust grain (Figure [2|. We per¬ 
form this simulation, and all other simulations, using the 
hybrid integrator of MERCURY suite of N-body integrators 
(Chambers 1999), which we have modified to include the ef¬ 
fect of radiation of dust grains, as given in equation[l] To test 
that it was implemented properly, we confirmed that our N- 
body code produces a K that is constant on long timescales 
with short timescale variations of size A K (equation [8| . 

We present here a particle from a simulation with a solar 
mass star, an earth mass planet on a circular orbit at 1 au, 
and a dust grain with /3 pr = 0.01, eo = 0.01, ao = 2.22au, 
i 0 — 0.0628. This particle is from the simulation labelled 
B in table [l] where we describe the total range of system 
parameters we consider in section |T2| The particle evolves 
initially under only PR drag, decreasing its semimajor axis 
as per equation [3] and its eccentricity as per equation [4] At 
~ 10 5 years, it encounters the 2:1 mean motion resonance, 
which does not capture the particle, but provides an eccen¬ 
tricity kick. Similar jumps are discernible when it crosses 
the 3:2, 4:3, and 5:4 mean motion resonances. 

At 1.5 x 10 5 yrs, the particle is captured in a 6:5 mean 


3.2 Suite of simulations: Setup 

The central star is characterised by one property, its mass 
(ra*). In all simulations, ra* = m©. 

We include one planet of mass ra p , with semimajor axis 
a p , with eccentricity e p , inclination i p , longitude of pericen- 
tre (tup), ascending node (H p ), and mean anomaly (M p ). In 
all cases, e p = 0, i p = 0. In all cases ti7 p , Q p , and the initial 
M p are chosen randomly and uniformly between 0 and 27r. 
Unless otherwise stated m p — m© and a p = lau. 

The dust grains are assumed to have zero mass, and 
the orbits are characterised by the same properties as the 
planet: semimajor axis a, with eccentricity e, inclination i, 
longitude of pericentre (zu), ascending node (Q), and mean 
anomaly (M). As with the planet, tu, Q, and the initial 
M are chosen randomly and uniformly between 0 and 27r. 
Unless otherwise stated, the dust grains begin with eccen¬ 
tricity] of 0.01, inclinatiorj^] of 2tt X 10 -2 , and semimajor 
axis chosen randomly and uniformly between 2.20 a p and 


2 This is the initial eccentricity of the dust grain, and neglects the 
consideration that high /3p r are created at higher eccentricities 
than their parent bodies (j |4.1| ). We incorporate that only when 
considering a full disk model (section [4| 

3 Comparisons to runs at i = 0.03147and i = 0.0157 show no 
significant differences. 
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Figure 2. Evolution of a dust grain in semimajor axis (top left), eccentricity (top right), resonance angle (bottom left), and resonance 
width (bottom right). The dust is affected by PR drag and radiation pressure, and perturbed by a planet with a p = lau and m p = m®. 
At 10 5 yrs, the dust grain crosses the 2:1 resonance and is perturbed to a higher e. Similar jumps are discernible when it crosses the 3:2, 
4:3, and 5:4 mean motion resonances. At 1. 5 X 10 5 yrs, the particle is captured in a 6:5 mean motion resonance with the planet. The 
eccentricity grows quickly in while in resonance, flattening out at a constant value, closely following equationThe dust grain escapes 
the resonance at 3.6 X 10 5 yrs. After escape PR drag brings the particle into the star at 3.8 X 10 5 yrs. 


2.25 a p . We choose this initial semimajor axis as we are in¬ 
terested in first and second order mean motion resonances, 
the most distant of which is the 3:1, which for a dust particle 
with /3 pr = 0 is located at 2.08 a p , and closer for /3 pr > 0 
(equation 10). Thus, a — 2.20 — 2.25 a p is effectively an 
arbitrarily large distance for all choices of /3pRp] 

We aim to be able to simulate disks over a wide pa¬ 
rameter space in /3pr, m p , a p , and eo- A complete sur¬ 
vey of the parameter space is computationally prohibitive. 
We choose a canonical simulation /3 pr = 0.01, mn p — m®, 
a p — lau, eo = 0.01 (simulation B), and explore these di¬ 
mensions by holding three of the parameters constant while 
varying the fourth on^] We consider /3 pr = 0.005...0.32, 
m p = l?n®...256ra®, a p = lau...l6au, eo = 0.01...0.64, in all 


cases spacing trial points by factors of two in the relevant 
parameter (table [l]) . Thus we perform a total of 25 simu¬ 
lations, each with 10 4 particles. Simulations are run with 
an 8 (0.01 / Pp r) (a/lau) 1 ' 5 day timestep, with data outputs 
every 10 4 (0.01 / Ppr) (a/lau) 1-5 days. Particles are removed 
from the simulation when they are more than 100 au, or less 
than 0.05 au, from the star. 

As all results are for ra* = M®, any fits to results may 
have unknown stellar mass dependence. Because our simu¬ 
lations generate a negligible amount of k ^ 2 captures, we 
only analyse the j + 1 : j captures in our simulations. 

Simulations continue until all particles have been re¬ 
moved. 


4 Because PR drag decreases e and a for particles with larger 
initial a, this isn’t strictly true. We can use equations [3] and to 
evolve dust grains at higher a to an e at 2.2au, so this choice is 
effective for our purposes. 

5 Non-zero initial e and i are more physical than exactly zero, and 
avoid possible degeneracies/singularities that might occur when 
they’re precisely zero; in particular, the resonant angle ip and its 


components (A P ,A, vo) are not well defined for zero eccentricity 
and inclination 
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Simulation 

eo 

/^PR 

m p 

a p 

A 

0.01 

0.005 

m® 

1 au 

B 

0.01 

0.01 

m® 

1 au 

C 

0.01 

0.02 

m® 

1 au 

D 

0.01 

0.04 

m® 

1 au 

E 

0.01 

0.08 

m® 

1 au 

F 

0.01 

0.16 

m® 

1 au 

G 

0.01 

0.32 

m® 

1 au 

H 

0.02 

0.01 

m® 

1 au 

I 

0.04 

0.01 

m® 

1 au 

J 

0.08 

0.01 

m® 

1 au 

K 

0.16 

0.01 

m® 

1 au 

L 

0.32 

0.01 

m® 

1 au 

M 

0.64 

0.01 

m® 

1 au 

N 

0.01 

0.01 

2m® 

1 au 

O 

0.01 

0.01 

4m® 

1 au 

P 

0.01 

0.01 

8m® 

1 au 

Q 

0.01 

0.01 

16m® 

1 au 

R 

0.01 

0.01 

32m® 

1 au 

S 

0.01 

0.01 

64m® 

1 au 

T 

0.01 

0.01 

128m® 

1 au 

U 

0.01 

0.01 

256m® 

1 au 

V 

0.01 

0.01 

1 m® 

2 au 

W 

0.01 

0.01 

1 m® 

4 au 

X 

0.01 

0.01 

1 m® 

8 au 

Y 

0.01 

0.01 

1 m® 

16 au 


Table 1 . A table of all the simulations we performed. Simula¬ 
tion B forms a centre, with other simulations varying by factors 
of 2 along axes of initial eccentricity, particle size (/?pr), planet 
mass, or planet semimajor axis. 


3.3 Suite of Simulations: Results 

3.3.1 Trapping probability 

To automatedly identify resonance capture, we look for dust 
grains that are continuously between aj.j+k + 2Aa max and 
aj:j+k — 2Aa max for at least twice the time predicted by di¬ 
viding the crossing distance 4Aa max by the instantaneous 
migration speed when the grain first encounters the reso¬ 
nance, given by equation [3] We employ the factor of 2 in the 
expected width to allow for deviations from our expression 
for the width owing to higher order terms. We also require 
the dust grain to be at that semimajor axis for at least twice 
as long to avoid falsely identifying dust grains as captured 
either when they are slowed by the resonant perturbation 
but not captured, or when oscillations (as equation |8| cause 
variations in a that result in dust grains staying slightly 
longer at the resonance than equation [3] predicts, or when 
close encounters with the planet cause a to evolve other than 
as predicted by equation [3] using the eccentricity of the dust 
grain when it first encounters the resonance. A spot check 
of 80 particles found that this approach correctly identified 
57 of 58 times when a particle crossing a j + 1 : j reso¬ 
nance was captured, and 359 of 363 times when a particle 
crossing a j + 1 : j resonance was not captured, considering 
only crossings up to the first reported capture. The trapping 
lifetime (or resonance escape time) is recorded as the total 
time the dust grain spends between aj.j+k + 2Aa max and 
aj..j+k — 2Aa max less the time predicted by integrating equa¬ 
tion [3] across that change in semimajor axis. 

We compare our observed capture probabilities to those 
calculated in §2.3| and plot a comparison of the cumulative 
capture probability as a particle crosses for different values 



5:4 10:9 15:14 


Resonance 


Figure 3. Cumulative capture probability for particles with 
eo 1 0.01, io = 27 r X 10 -2 , a p = lau, m p = m®, for various 

/3pr (simulations A-G). Comparing the outcomes of our simula¬ 
tions (red solid lines) with the calculations in |2.3| (blue dotted 
lines). The fit works well overall, although the critical j cutoff is 
perhaps smoother in the N-body simulations than in our Hamil¬ 
tonian approach. 



Figure 4. Cumulative capture probability for particles with 
/3 pr = 0.01, io = 27 r X 10 -2 , a p = lau, m p = m®, for var¬ 
ious eo (simulations B, H-M), comparing the outcomes of our 
simulations (red solid lines) with the calculations in j |2.3| for dif¬ 
ferent initial eccentricities (dotted blue lines). Each initial eccen¬ 
tricity is offset vertically by 0.75 for clarity. The eo = 0.32 and 
eo = 0.64 case both have < 1% capture in both N-body simula¬ 
tions and the predictions of §2.3| and are not plotted here. At low 
eccentricity the agreement is very good, and it remains reasonably 
good at all eccentricities. 


of /3 pr (figure p - simulations A-G) and for different val¬ 
ues of eo (figure 31- simulations B, H-M). This comparison 
considers only the first capture of particles in the N-body 
simulation; any subsequent captures are neglected as eccen¬ 
tricity growth while in resonance (£2.5) makes such com¬ 
parisons unsuitable for evaluating particles with the initial 
conditions of the experiment. This comparison accounts for 
the eccentricity evolution as the semimajor axis decreases, 
and the eccentricity kicks due to resonance crossings (figure 
[2|; the quoted eccentricities are the initial values. The calcu¬ 
lated capture probabilities match the model well, although 
slight differences exist. 

The 2:1 resonance presents a particular challenge as 
asymmetric libration results in capture into both the 2:1(1) 
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Figure 5. The critical j above which capture cannot occur. Lower 
bounds come from simulations A, B, and C, where the N-body 
simulations capture ~ 100% of particles, so capture must be able 
to occur until the j where the cumulative capture probability is 
~ 100%. The upper bound is from simulation F, where the lack 
of any captures constrains the maximum j to be less than the 
lowest j for which capture is predicted in the Hamiltonian model. 
The N-body estimates are from simulations D and E, with the 
critical j taken from where the Hamiltonian model and N-body 
simulations predictions for capture probability diverge. The green 
points show the j given by equation |27| with A = 2.3 chosen to 
bring the Hamiltonian model into agreement with the N-body re¬ 
sults. The solid pink shows where the dust grain would orbit at 
the same semi-major axis as the planet; here the impact parame¬ 
ter goes to zero, and hence the calculated kick becomes infinitely 
large. 


and 2:l(u) resonance, with differing probabilities (Message 


1958 Chiang &; Jordan|2002 Murray-Clay &; Chiang|2005 ). 


We emulate the approach of Wyatt (2003) and calibrate the 


relative capture probabilities against our simulations in the 
form 


P 2 -.i (/) = 0.5 - a6>V, (33) 

where 0 = (a/lau/Myr) /(a/lau) (Mg/m*). We per- 
formed limited extra simulations of capture where capture 
into both the (u) and (1) resonance occurred, with 1000 
particles; sampling planet mass every 25ra® from 150 — 
400m® (with (3 — 0.01) and migration rate with factors of 
2 in (3 from j3 — 0.0025 to (3 = 0.02 with m p — 256m®. Our 
best fit parameters were a = 0.01, b = 0.25, and c = —0.4, 
which we employ in ^4]to predict the relative rates of capture 
into the 2:1(1) and 2:l(u) resonances. 


3.3.2 Maximum j cutoff 


Our derivation of the highest j at which capture can occur 
has a numerical coefficient A that was not specified in §2.4| 
Comparison of our Hamiltonian capture model with our N- 
body simulations favours a coefficient of A « 2.3 to bring 
the two into agreement (figure [5]). This is slightly higher 
than values previously reported (see section 2.4). 


3.3.3 Libration width at capture 

The model discussed in section |2.3| makes a prediction for 
what the distribution of libration widths should be at cap¬ 
ture (Spo). In figure [6] we compare that prediction to our 


Figure 6. Cumulative distribution of the resonant angle libration 
width at capture (<fy?o)- Each resonance is offset by 1 vertically 
from one another for clarity. Width is measured using the max¬ 
imum and minimum libration angles (cp) while 0.04 < e < 0.06. 
Thick lines are simulations, thin lines are predictions from sec¬ 
tion |2.3| (for the most part, the two lines overlap too much to 
be distinguished). These are all taken from the default case of 
(3 pr = 0.01, eo = 0.01, rrip = m®, a p = lau (Simulation B). The 
overall agreement is excellent, except in the 4 : 3 case, where the 
simulation had only 55 captures. 


simulation results from simulation B, the canonical case, us¬ 
ing the maximum and minimum libration angles (<p) while 
0.04 < e < 0.06. This is chosen to be confident the particle 
is in resonance, but has not experienced significant libration 
width growth. The agreement is generally excellent, except 
in the lowest j case, where the predicted libration widths are 
slightly too large. We will see that it is important this dis¬ 
tribution is well characterised, as the initial libration width 
sets the length of time the particle remains in resonance 
after capture. 


3.3.4 Eccentricity evolution 

In section |275l we presented an expression for the eccentricity 
evolution of a particle caught in resonance (equation [M|. 
Here, we evaluate how closely that expression matches what 
we find in N-body simulations. 

In figure [7| we compare the evolution of eccentricity 
while in resonance for particles from simulation B to equa- 
tion[3T] plotting all particles caught in 5:4, 6:5, 7:6, 8:7, 9:8, 
and 10:9 mean motion resonances for at least 5 x 10 4 years. 
The evolution timescale is in excellent agreement, and the 
magnitude is in reasonable agreement, with only a ~ 4% de¬ 
viation of the simulation from the analytic solution. The 
analytic solution assumes particles are captured at zero ec¬ 
centricity, but as the growth is an exponential decay away 
from 0, this difference is not significant. 


3.3.5 Centre of Libration 

As seen in figure [2] after capture, the eccentricity rises, and 
the predicted centre of libration moves quickly towards a 
value set by the final eccentricity. The centre of libration 
found in simulation lags the predicted value, however, as 
seen in the bottom left panel of figure [2] The same general 
behaviour of an initial offset with a linear relaxation towards 
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Initial libration width 
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Figure 7. A comparison between the expected evolution of the 
eccentricity e prec iicted predicted from equation [31] (solid lines), 
and that found in the simulations e s i m (dots). All particles plot¬ 
ted here are from the canonical case of /3 pr = 0.01, eo = 0.01, 
m p = ra®, and a p = lau (simulation B). Plotted are 1000 ran¬ 
dom particles caught in 5:4 (black), 6:5 (red), 7:6 (green), 8:7 
(blue), and 9:8 (magenta) mean motion resonances for at least 
5 X 10 4 years. 


the equilibrium value is also found for the other particles in 
the simulations. Thus, we set out to characterise this be¬ 
haviour with a phenomenological model of the form 

<Po irn = (l - — ) . (34) 

\ Trelax J 


Here the maximum offset in the centre of libration, Ac^o and 
the relaxation time, r, are parameters that we will fit phe¬ 
nomenologically, and t is the time after resonance capture. 

We begin with a fit of the maximum offset (A<^o) be¬ 
tween the equilibrium and actual centre of libration, from 
simulation B, plotted against the initial libration width, for 
several resonances (figure [8]). Here, the maximum offset is 
found to depend only on two parameters, the j of the reso¬ 
nance, and the initial libration width of the captured parti¬ 
cle, dp o (j ]3.3.3 ). The best fit is of the form 


Ap 0 = Ci cos (5ipo/C2). 


(35) 


Fitting the same for all simulations, we find 

/ \ —0.864 _ n 49 o 

Ct = 4475/#rr°- 81 (^j) (59 " ' (36) 

and 

C 2 = 163.9 - 1.76j - 1.4 (+ 0.73 (WV (37) 
V ?7i® J Vlau/ 


We neglect the /3 pr dependence of C 2 as the best fit coeffi¬ 
cient of —1.93 d= 1.76 is consistent with zero. 

Plotting the maximum offset against the relaxation 
time, we find the two have an approximately functional cor¬ 
respondence (figure [9]). We fit the maximum offset between 
the predicted centre of libration and that found in simula¬ 
tions against time it takes for the centre of libration to to 
relax the equilibrium value of equation [29] The best fit value 
is 


Figure 8. Maximum offset between the predicted centre of li¬ 
bration and that found in simulations, compared to the initial 
libration width, for particles from simulation B. The solid lines 
are a fit to all simulations, given by equations |35| [36| and |37| 



200000 300000 400000 

Relaxation time (years) 


Figure 9. The time it takes for the libration centre to relax to the 
equilibrium value of equation [29] compared to the initial libration 
width, for particles from simulation B. The solid lines are a joint 
fit to all the simulations of the form Apo = C 3 (l — e~ Trelax / CA ^j. 


with 

C 3 = 7605j“ 103 /3pRmp°- 94 a“ 0 - 45 (39) 

and 

Ch = 13949j p PR m a p . (40) 

An example of this fit for simulation B is plotted in figure 
[9] Here, of course, the combination of equations [35] and [38] 
allows one to express the relaxation time as a function of 
the initial libration width in the form 

Trelax = C 4 In (f - A COS A , (41) 

and thus for particles with the same system parameters 
caught in the same resonance, their evolution is dictated 
entirely by their initial libration width. 


3.3.6 Libration width evolution 


Trelax 


= Chin 



Ap 0 \ 
C 3 J 


(38) 


The resonance width (Sip) was found to grow exponentially 
with time by Sidlichovsky & Nesvorny (1994), but they did 
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Figure 10. The e-folding time of the resonance width for 
dust grains in all simulations. We fit the resonance width 
as exponential growth between 160(0.01//3pr) (a p / lau) 2 kyrs 
and 320(0.01/^pr) (a p /lau) 2 kyrs for those particles which are 
trapped in resonance for at least 480 (0.01//3pr) (a p / lau) 2 kyrs 
(all simulations not plotted have zero particles meeting this cri¬ 
terion). Points are values from the N-body simulations, the solid 
lines are equation |43| Here the error bars cover the smallest to 
largest values, with the point centred on the median value. Rare 
catastrophic failures of the automatic fitting produce extremely 
large values, so we employ the median value, rather than the mean 
value. 



0 0.2 0.4 0.6 0.8 1 


Time (resonance escape times) 

Figure 11. Closest approach between the planet and the dust 
particle in simulation B. Dust grains are binned with all grains 
trapped in the same resonance in resonance escape time bins with 
width 10 4 yrs. Closest approach is sampled across 1000 years 
for each bin. Once the particles come within ~ 0.04 au, they 
are ejected from the resonance. We measure this escape radius 
across all simulations using the second last time bin, where 10% 
of measurements are below the escape radius. Only those particles 
caught for more than 10 4 years are included. 


not quantify this behaviour. We also observed exponential 
growth with time (Figure [5]). This growth is of the form 

Sip — (42) 


where 5ipo is characterised in §3.3.3| To fit the reso¬ 
nance width of a particle in resonance, we bin all out¬ 
puts of ip within a 4000 (0.01/^pr) (a/1 au) 1 ' 5 year pe¬ 
riod (as in £3.3.5), taking the largest difference between 
two values as the resonance width at that time. For those 
particles that persist in resonance for more than 4.8 x 
10 5 (0.01/^pr) (a/1 au) 1 ' 5 years, we fit an exponential curve 
to the growth between 1.6 x 10 5 (0.01//3pr) (a/1 au) 1 ' 5 and 
3.2 x 10 5 (0.01//3pr) (a/1 au) 1 ' 5 years to assign a to that 
particle’s evolution. The of different particles with the 
same /3pp caught in the same resonance cluster tightly (Fig¬ 
ure 


10 ). 


By fitting all of the values derived this way for difference 
resonances and simulations, we find 


T ' p w L14 x 105 O Crr) (£) 2years ’ (43) 

a fit of which can be seen in plot |10| Fixing the constants 
of equation |43[ varying one while keeping the rest fixed, a 
non-linear least-squares Marquardt-Levenberge fit, weight¬ 
ing each point by y/n where n is the number of particle- 
histories represented by a point, gives the normalisation as 
1.14 =t 0.01 x 10 5 , the exponent of the (j + 1) /j term as 
2.02db0.01, the exponent of the /3pp term is —0.94±0.05, the 
exponent of the planet mass term would be —0.006 d= 0.002 
if we included a power law term dependent on the planet’s 
mass, and the exponent of the semimajor axis is 1.99 =L 0.02. 

Note that the scalings in equation ^3] are very similar 
to the PR drag lifetime in equation [5] They are related to 


one another as 


3 + 1 


(1 — /3pr) 3 Tpr 


(44) 


and here the j and /3 pr terms will be close to one, and thus 
the e-folding time will always be larger than the PR drag 
time by a factor of a few. Because any resonant object must 
have 5ip < 7r, unless dust grains are captures at very small 
libration widths (which is not found by |Mustill &; Wyatt 
(2011) for any parameter space they investigate), they can 
only undergo a few e-folds of their libration width before 
they escape the resonance. Therefore, resonantly captured 
rings should never exceed the background disk in brightness 
by more than a factor of O (10). 


3.3.7 Escape from resonance 

At early times, dust grains undergo close encounters with 
the planet as they are pushed into the resonance. The res¬ 
onance then protects the dust grain from close encounters, 
but as the eccentricity and libration width increase, encoun¬ 
ters become increasingly closer, until the particle is lost from 
the resonance. To analyse if and when close encounters cause 
a dust grain to be lost from resonance, we bin together all 
dust grains from the same simulation which are caught in 
the same resonance and which are trapped in that reso¬ 
nance for the same amount of time (rounded to the near¬ 
est 10 4 (a/1 au) 2 (m p /m®) 0 ' 5 years). We define the closest 
approach for each bin of dust grains by the smallest planet- 
dust grain separation output by the simulati on over each 
10 3 (a/1 au) 2 (m p /m®) 0 ' 5 year interval (Figure[ll|. Because 
the planet-dust grain closest approaches decrease during the 
evolution, and escape occurs for all particles at a fixed en¬ 
counter distance, we conclude escape from resonance occurs 
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Figure 13. Lifetime of particle in resonance plotted against the 
time it takes the centre of libration to relax to the equilibrium 
value. Points are from simulation B. A fit to all simulations of the 
form T r es = t re lax + C is plotted for the appropriate resonances. 


due to a close encounter with the planet, which provides an 
impulse sufficient to move a particle out of resonance. 

To measure this distance, we consider the 10 th per¬ 
centile closest approach in the second to last time bin in 
the separation (to allow for some error in the automated fit¬ 
ting, sampling time, and identification of when escape from 
resonance occurred). Because eccentricity evolves to a fixed 
value (§3.3.4) we do not expect eccentricity at capture to 
affect the resonance escape distance, and our simulations 
are compatible with no e dependence (figure 12). Addition¬ 
ally, no /3 pr, nor j dependence is observed. The escape ra¬ 
dius is observed to depend most strongly on the planet mass 
and semimajor axis. A non-linear least-squares Marquardt- 
Levenberge fit of the escape radius in simulations B and N-Y 
gives the escape radius as 

0.616±0.056 


R e « 0.0361“ 


of the planet. 


V 


11 ip 

m@) 


( a P \ 

V lau / 


0.931±0.012 


(45) 


3.3.8 Resonance lifetimes 

Comparing the relaxation time to the total time spent in 
resonance shows that the two show a functional relationship, 
which can be fitted as a linear relationship with a constant 
offset that depends on the resonance and other properties of 
the simulation. We perform a joint fit across all simulations 
of 


Tres — Trelax 


+ c 5 , 


(46) 


Figure 14. Time in resonance versus the initial libration width 
in simulation B (points). The resonant lifetime is fixed for a given 
capture by the initial libration width and j, as the libration width 
grows predictably in that case (equation |43| ), until the particle 
has a close encounter with the planet (equation |45|), causing it 
to be ejected from resonance. This results in a predictable res¬ 
onance lifetime (the lines are equation |48| for the corresponding 
resonances). 


time in resonance as r res = Ca In (1 — C\ cos (Scpo /C^/C^T 
C5. In practice, r res changes very rapidly at small ipo, and 
the uncertainties accumulated in bootstrapping make the 
expression a poor fit. Using the general expression, we refit 
with 


Tres = Ca In ^1 - cos + Cc- (48) 

In equation [48] the best fit values across all simulations are 

C.--1.27 x (a) (£)',«. 

(49) 

C» = 66-4.4, -3.2 (=£)- 0.7 (i), (50) 

and 

/ \ 0.01 \ 9 n 

Co = 3050,-° W (i) (£)• (51) 

years. A comparison of this fit to the results of simulation 
B is plotted in figure [14] 

4 APPLICATION: MODELLING A WHOLE 
DISK 


and get a best fit of 


C 5 = 377.5 j 66 p 


•-1.33/3-1.54 0.21 2.2 


'PR 






(47) 

An example of this is plotted in figure [I3| for simulation B. 
Knowing how dust grains escape from resonance (equa- 
45), and how they evolve in resonance (equations |31| and 


tion 


43), we can, in principle, predict the resonance lifetime of a 
trapped dust grain. That analysis does not, however, lend 
itself to a neat analytic expression. Instead, we derive a life¬ 
time by combining equations [35] [38] and [46] which gives the 


We are now in a position to predict the dynamical evolution 
of a population of dust grains as they evolve past a planet. 
Disk images are generated by the following step-by-step ap¬ 
proach: 
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Figure 12. Escape distance vs. planet semimajor axis in simulations B, V, W, X, Y (top left), escape radius vs planet mass, from 
simulations B, N-U (top right), escape radius vs /3 pr, from simulations A — C (bottom left), and escape radius vs initial eccentricity, 
featuring simulations B, and H through M (bottom right). Error bars are the 10 th , 20 th , and 30 th percentile fits to the second last time 
column in plots like figure [TT] for the relevant simulations. The fit (blue dotted line) is a joint fit to mass and semimajor axis simulations 
(B, N-Y). 


1. Initial conditions are specified by assigning a mass to 
the star, mass and semimajor axis to the planet, and semi¬ 
major axis, eccentricity, and inclination to the parent body 
of the dust grains, and /3 pr to the dust grains. To consider 
a population of parent bodies (as we will in ^5|, the parent 
body for each grain is chosen randomly from the population. 

2 . Dust grains are placed on orbits that correspond to 
their parent body orbits (modified for their /3pr, as outlined 
in §4.1) To consider a population of grains with different 
/3 pr is somewhat complicated, as the number of grains may 
be a steep function of /3pr. This can result in under-sampling 
or oversampling, and we create images for single /3 pr values, 
and weight and add the images afterwards to model a grain 
population. 

3. Dust grains evolve in ad and ea due to PR drag as 
calculated by Wyatt k Whipple ( 195Q| (j |2.1| ). 

4. Whenever a particle encounters a first order mean mo¬ 


by the model of Mustill k Wyatt (2011), updated in £2.3 


tion resonance with the planet (§ 2.2), it may be captured 


into resonance with probability given in £2.3 (step 5.A) or 
not (step 5.B). 

5.A.i If a particle is captured into resonance, the ini¬ 
tial libration width is taken from the distribution given 


The initial centre of libration is offset from the equilib¬ 
rium value of equation [29] by an amount dependent on 
the resonant width at capture, as given by equation [35] 
5.A.ii While trapped in resonance, particles evolve in 
eccentricity as per equation |3 1 1 and libration width as per 
equation |42| The offset from the equilibrium centre of li¬ 
bration relaxes linearly to the value given by equation |29| 
in a relaxation time given by equation |38| 

5.A.iii When the particle comes within a distance 
given by equation ^5] of the planet, it is removed from 
resonance, and begins evolving again under PR drag (i.e., 
return to step 3.). 

5.B If a particle is not captured into resonance, an ec¬ 
centricity jump is applied as found in |Mustill k Wyatt| 
(J011) (with the modifications found in the Appendix), 
and the particle continues evolving under PR drag (i.e., 
return to step 3.). 

6. When the dust particle encounters the star it is re¬ 
moved from the simulation. 

We plot comparisons of the N-body disks of table [l] 
to images produced using this semi-analytic prescription 
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in figures [15] and |16| N-body disk images were produced 
by sampling the location of the N-body particles every 

_ i 

10 3 /^pr a p 5 days, where a p is measured in au. To produce 
an image, at each time we rotated the dust grain around the 
star so the planet is located at (1,0). The semi-analytic disks 
are very good matches for the N-body disk, with increasing 
discrepancies for higher mass planets. In these cases, cap¬ 
tures into the 2:1 resonance are not reproduced as well, as 
the grains move between the eccentricity and inclination res¬ 
onances, which we do not model. Additionally, while there 
is reasonable agreement with our model in the parameter 
space we probe with iV-body simulations, as the calibration 
of capture into the (u) and (Z) 2 : 1 resonances in equation 
|33| is calibrated over that range, it is unclear how far it can 
be extrapolated outside this range. 


4.1 Collisional production 


A real disk will contain particles with a variety of initial 
eccentricities, and a variety of different Apr- These two dis¬ 
tributions must be modelled to produce a realistic disk. Of 
note, the initial dust properties will differ from the parent 
bodies, owing to radiation pressure, so the dust grains will 
begin with a distribution of semimajor axes and eccentrici¬ 
ties that differ from those of the parent bodies and are func¬ 
tions of Apr (i.e., ad (Apr) and ed (/3pr)). If a collision oc¬ 
curs at r, neglecting the velocity at which dust is ejected 
from the collision, the dust orbital parameters relate to the 
parent body parameters (a b and e b ) as: 


dd 


£d 


( _1_ 

\ (1 — Apr) CLb 


1 — /3pR 


2/3pr \ 1 

(1 - /3pr )r) 


^1 + 2 Apr 


1 - el a b 
(1 — /^pr) 2 r 


ApR 

(1 - Apr) 


1-eg 

(1 - Apr) 2 


2 


(52) 

(53) 


where the approximate expressions assume the parent bod¬ 
ies’ semimajor axis a b ~ r and the parent bodies’ eccen¬ 
tricity e b ~ 0. This last result motivates the choice of 
eo = Apr as a default case. Given the parent body pop¬ 
ulation and the distribution of r for generative collisions, 
the initial orbits of particles as a function of Apr can be 
predicted (Wyatt et al. 2010) (e.g., Figure [Tt| shows the ex¬ 
pected dust starting properties from collisions of asteroids 
in main belt in the Solar System.). 


5 EARTH’S RESONANT RING 


In addition to the general ring structure produced by reso¬ 
nantly trapped dust ( Kelsall et al.f l998), a smaller pattern 
exists in which a ‘blob’ of dust trails the Earth in its orbit 
(Dermott et al. 1994; Reach et al.| [l995). Measurements of 
the North Ecliptic Pole sky brightness by the Spitzer Space 
Telescope, which is on an orbit that causes it to slowly drift 
away backwards along the Earth’s orbit, have measured the 


azimuthal structure of this blob in 8/xm emission (Reach 


2010 ). 



Eccentricity 

Figure IT. Initial eccentricity of dust particles compared to a 
parent population. The parent population’s eccentricity distribu¬ 
tion is that of the first ten thousand numbered asteroids. The 
distance from the sun at which the collision takes place (r) is 
chosen uniformly between periapse and apoapse; a real popula¬ 
tion would be expected to be biased towards periapse (|Wyatt] 
|et al.p blP]). The two distributions are roughly the same at low 
/3 pr, and diverge when /3 pr ~ e b . 


We apply our model to the structure of the trailing blob. 
For the orbits of the parent bodies, we consider the first 
10000 asteroids with orbital elements taken from the Minor 
Planet Centre’s Orbit (MPCORB) Database. We sample the 
size distribution at Apr = 0.45, 0.4, 0.35, 0.3, 0.25, 0.18, 
0.13, 0.088, 0.07, 0.063, 0.044, 0.031, 0.022, 0.016, 0.013, 
0.011, and 0.0055. The relative components are calculated 
separately, and added together to produce the spatial den¬ 
sity of small grains. This requires a size-number distribu¬ 
tion of grains, and here we employ the result of |Love fe| 


Brownlee (1993), who counted the craters from dust im¬ 


pacts on the space-facing side of the Long Duration Ex¬ 
posure Facility, and converted that to a size distribution 
of dust using experimental data on the impacts of dust 
grains onto aluminium targets and assumptions about the 
impact velocity. We equate particle sizes with Apr per equa- 
tion [2] assuming p = 2.5 g cm -3 . Here a concern exists; 


Love & Brownlee (1993|) only measured crater sizes between 


20/im and 1400/im, corresponding to f3 < 0-07 or s > 3/xm. 
At smaller sizes, their best fit size distribution employs a 
polynomial which turns up strongly in a way that may not 
be physical; we instead affix a power-law distribution be¬ 
tween s — 3/xm and the blow-out size (figure [l8| . 

With this we can then predict the light curve Spitzer 
would see on its actual path using our model disk. We ob¬ 
tained Spitzer’s trajectory from JPL’s HORIZONS system 
(Giorgini et al. j! 1996). Luminosities are calculated by treat¬ 
ing the dust grains as spherical black bodies 3 ! and integrat- 


6 With most of the emitting area in large grains (figure [l8] ), this 
should be fairly reasonable. For the smallest grains, this may get 
the temperature of the smallest grains wrong by less than a fac- 
tor of 2 (|Backman Sz Pa resce 1993]) . Using the grain models of 
|Augereau et al.| ( |1999| ), as implemented by |Wyatt fc Dent| ( |2002| ), 
we estimate the total flux as a function of size might be uncer¬ 
tain by a factor of a few at the smallest sizes, which would need 
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Figure 15. Simulations A - G (variations of /3 pr) left (1st column is semi-analytic approach, 2nd column is N-body approach), and 

Simulations B, H-M (variations of eo) right (3rd column is semi-analytic approach, 4th column is N-body approach). In all plots, 10000 
... .... _ —n * i c . . . . __ _ .. . „ „ „ m 
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Figure 16. Simulations B,N-S (variations of m p ) left, and Simulations B, V-Y (variations of a p ) right (and T-U, variations of m p ), 
bottom right. In all plots, 10000 particles are each plotted once every 10 3 X /3p^' 5 a 1-5 days, in a grid of 400 by 400 cells covering from 
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Figure 18. Best fit size distribution of dust grains near the Earth, 
as measured by |Love &; Brown lee (1993) (thin red line). The thick 
red line segment is the size distribution over which their measure¬ 
ments were made. As the upturn below ~ 3p,m is not attested 
from the data, we use an alternate size distribution at small sizes, 
the blue line (of intermediate thickness). 



Figure 19. Zodiacal cloud luminosity expected to be observed by 
Spitzer from our model of asteroidally produced dust (thick red 
line), compared to measurements by |Reach (2010| (blue points). 
A constant background with an annual sinusoidal modulation has 
been subtracted off; the remaining flux is normalised such that 
the component subtracted off has a mean amplitude of 1. The 
model clearly does not match the data. 


ing the expected emission along the line of sight of Spitzer. 
To compare with measurements of the trailing blob, we ap¬ 
ply the same approach as Reach (2010); a constant ampli¬ 
tude signal, with a sinusoidal modulation is fit to the data, 
and subtracted off. We plot the residuals in figure [19] The 
predicted amplitude of the signal is then much stronger than 
the observations. 

Detailed inspection of the contributions from different 
/3 pr finds that the disagreement comes mostly from the 
larger grains, which are trapped at low j (figure [20] and re¬ 
call figure [T5| . This produces a large dip at the Earth, that 


to be included in a more detailed model. But as they contribute 
relatively little flux, we do not attempt to model that here. 


Figure 20. Zodiacal cloud luminosity expected to be observed 
by Spitzer from our model of asteroidally produced dust of three 
different mono-size distributions (Jpr = 0.4,0.063, and 0.0055), 
compared to measurements by |Reach| |2010) (blue points). A con¬ 
stant background has been subtracted off; the remaining flux is 
normalised such that the component subtracted off has an ampli¬ 
tude of 1. From this, we can infer that the addition of more small 
grains, or grains that are not trapped into resonance for other 
reasons, could bring the model in figure [T9| into alignment with 
the data. 


rises to a flat level at later times. As such, the model can be 
brought into rough agreement with the data by the addition 
of a smooth disk component. We plot two such examples in 
figure |2lj In one example, we simply use asteroidal grains, 
but do not correct the Love & Brownlee (1993) polynomial 
fit to the size distribution at small sizes where the size distri¬ 
bution extrapolation is suspect. In this case, the high migra¬ 
tion rate of the high /3 pr grains prevents their capture (recall 
figure [l]- such grains have high dB/dt). In another, as it has 
been suggested the zodiacal light may come mostly from 
cometary grains ( Nesvorny et al.||2010| , we assume 80% of 
the zodiacal light comes from comet ary grains, which we 
model using Halley’s comet as the parent body orbit, with 
all grains released at pericentre. These grains are not cap¬ 
tured due to their high eccentricity (again recall figure [l]- 
such grains have high Jo). Because of the high initial eccen¬ 
tricity, grains with /3 pr >0.01 are produced with e > 1.0, 
and so only grains with /3 pr ^ 0.01 are used. However, a 
more complicated population of parent bodies and ejection 
velocities would be needed for the cometary component to 
reproduce the Love & Brownlee (1993) size distribution, or 
a model in which asteroids and comets contribute differently 
to the size distribution at different sizes. Thus for now, this 
should be considered only a demonstration of principle, and 
we defer a more detailed model to a later paper. 

To aid in the visualisation of the structure of the ring 
and Spitzer’s path through it, we plot the azimuthal asym¬ 
metries of an 80% cometary grains and 20% asteroidal grains 
disk in figure |22[ along with the path of the Spitzer Space 
Telescope. 
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Figure 21. Zodiacal cloud luminosity expected to be observed by 
Spitzer using the uncorrected dust size distribution fit of | Love fc] 
|Brownlee| ( |1993| > assuming only asteroidally produced dust (solid 
red line), and using the corrected dust size distribution fit of | Love | 
| fc Brownlee| JL993), but assuming 80% of dust is produced by 
comets (dotted green line), compared to measurements by | Reach] 
(2010| (blue points). The cometary model assumes all dust starts 
with the orbit of Halley’s comet, and all dust is released at peri- 
centre. A constant background has been subtracted off; the re¬ 
maining flux is normalised such that the component subtracted 
off has an amplitude of 1. The model produces rough agreement 
with the data. 
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Figure 22. Disk flux relative to the azimuthal average at the 
same radius in our 80% cometary, 20% asteroidal grains model. 
The positions of the Sun and the Earth are plotted with ★ and 
® respectively. The path of Spitzer is plotted, showing its move¬ 
ment into Earth’s trailing blob. This plot also reveals that the 
trailing/leading asymmetry is the main asymmetry in the neigh¬ 
bourhood of the Earth, although a large asymmetry is just visible 
at ~ 2.2 au, the apocentre of low /3 pr particles in the 2:1 reso¬ 
nance. 



6 DISCUSSION 


Exozodiacal dust is a possible source of obscuration for fu¬ 
ture missions attempting the direct imaging of terrestrial 
exoplanets ( Beichman, Woolf &; Lindensmith||1999 Defrere 


et al. 2010). However, capture of dust into resonance can 


allow for the detection of a planet’s presence, even where di¬ 
rect imaging is impossible (Kuchner & Holman 2003). Quick 
generation of disk images is likely to be key to determining 
planetary properties from the appearance of a resonantly 
captured ring. In this work, we present a semi-analytic ap¬ 
proach for calculating the evolution of dust particles moving 
past a planet due to PR drag. This allows dust evolution to 
be calculated on the PR time-scale, rather than the orbital 
timescale, which allows orders of magnitude faster calcu¬ 
lations than traditional N-body methods. We compare this 
approach to results generated from N-body simulations, and 
find them to be an excellent match for terrestrial and super- 
Earth type planets (m < 10M©), which radial velocity re¬ 
sults ( Mayor et al.||2011| and transit results (Howard et al. 
2010) find to be the most prevalent types of planets. Larger 


planets capture many dust grains into the chaotic 2:1 reso¬ 
nance, for which our approach produces a reasonable match. 

We apply this model to dust produced in the Asteroid 
belt that is resonantly trapped as it passes the Earth. Using 
the measured size distribution of dust grains near the Earth, 
the model does not match the observations, requiring an ad¬ 
ditional smooth component to the Zodiacal cloud to bring 
the model into agreement with the data. This smooth com¬ 
ponent can be produced by comets, for instance, which re¬ 
lease dust grains with much higher initial eccentricities. Al¬ 
though other models can be considered, this could be taken 
as support for other lines of evidence that suggest the Zo¬ 


diacal cloud is mostly produced by comets, not asteroids 
( Nesvorny et al.||2010 ). 

Application of this model is likely to find great utility 
in exozodiacal systems, where the decrease in computation 
time can make it more feasible to address the large param¬ 
eter space uncertainly associated with the unknown orbits 
of parent bodies, and the unknown orbits and masses of 
planet. This improvement in computation time for study¬ 
ing the dynamical aspects of the evolution will also make it 
feasible to consider more realistic prescriptions for the col- 
lisional evolution of the debris as it evolves from the source 
location. While it was possible to ignore this for our study of 
the Solar system because the size-mass distribution of dust 
at the Earth has been measured directly, this would have 
to be calculated for the simulations for extrasolar systems. 
We intend to address this question in an upcoming work, 
and explore the application of this model to measuring disk 
structure with repeated LBTI observations of zodiacal light, 
or possible direct imaging by future missions. 
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APPENDIX A: MODIFICATIONS TO MUSTILL 
& WYATT 2011 FOR RADIATION FORCES 


To adapt MW11 for radiation pressure, we non- 
dimensionalise as in MW11 to put time in units of the 
planet’s mean motion and length in units of the planet’s 
semi-major axis. 

The Poincare variables for the particle will take the form 

A = \[B~a (Al) 

T = A (l - s/\ - e 2 j ~ Ae 2 /2, (A2) 

for some value of B. The Keplerian part of the Hamiltonian 
will take the form 

«^K e P = “^2 ( A3 ) 


for some A. These will not be the same as in the pr = 0 case 
because the radiation pressure changes slightly the location 
of the resonances, and hence the coefficients in MW11, but 
A as a function of /3 pr must tend to 1 as /3 pr 0. We 
can exploit the fact that the equations of motion must be 
canonical, so that the mean motion 


dJtfke p 

but we know that also 


A_ 

A 3 ’ 


(A4) 


n = a/(1 - f3p R )/a 3 (A5) 


and so 


A 2 B 3 = 1 — ^pr. 


(A6) 
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at the resonance location 


ao = (l-/5 PR ) 1/3 (t±i) 

(1-/3pr) 1/6 


— R 1 / 2 ( 


1/3 


Ao = B 


Expand about the resonance location with A = Ao + / 

A , A 3A 2 

^ Kep “ 2A 0 + Ap 2 aJ 

and drop the first term as it is constant. Then 
jB~ 3/2 (1 - 

k 4/3 


^Kep 




■2/3 


j +: 


— Ao 

where 


(A7) 

and 

a = ^(1 — /#p R ) 1/3 ao, 



-/3p R ) 1/2 aU0' + l)-3(l-/3 PR ) 2/3 

(A8) 

This latter determines the migration rate. 

We now make the usual transformations 

(A9) 

o' = 

f = 

7T — 6 , 

XJ U 


t! = 

Yt , 


b' = 

Zb , 


M” = 

-VEJT, 


The resonant term will have the form 
J% es = rr 1/2 cos [(j + 1)A - jApi - zu] 


(A10) 


(AH) 


for some r, with A and zu the mean longitudes and longi¬ 
tudes of periastron. We need to set the value of r. Lagrange’s 


Equations for a first order resonance give (Murray & Der- 
mott|1999 ) 


tt=~ + 1 


For symplecticity we require 

dT _ dJ4? res _ ^pl/2 
dt dO 
But from the definition of T we also have 


rT i/2 sin <9. 


(A12) 


(A13) 


dT de 

— = Ae —, 
dt dt 


(A14) 


and substituting in de/dt and equating coefficients gives 

V/3lffi) 


r = -2 1 ^B 7 /\f 31 a- 5/4 A- 1 +2- 1/2 B 7 ^al /4 tiA-\ (A15) 


which is Eq A7 of MW11 with some extra factors of A and 
B. 

Here we impose A = 1, since it is always possible to scale 
a Hamiltonian so long as time is also rescaled to maintain 
symplecticity, and here we have the correct time scaling via 
the mean motion. Fixing A — l gives B = (1 — /3pr)~ 1 ^ 3 
and hence the full Hamiltonian 


= 


(1 — /?Pr) 2 j _ 3(1 — /3pr) 3 ^ 


3/2 


2a 2 0 


2 2 m/3i(1 — /3pr) I 2 


— 2 2 (1 — /3pr) 12 CLq IjL 


x r 1/2 cos[(j + 1)A - jApi - zu] 


(A16) 


which is the revised Eq A14 from MW11. Next we make the 
usual transformations 


r = Ji, 

0 \ = (3 + 1)A - jApi - 07, 

The Hamiltonian becomes 


I = jJi + J 2 j 

6>2 = A. 

Jf 7 = —aJi + bJ 2 — rj\/ 2 cos 0\ 


(A17) 

(A18) 

(A19) 


(A20) 


(A22) 
(A23) 
(A24) 
(A25) 
(A26) 

where comparison of coefficents and enforcement of the 
canonical condition on the equations of motion requires 


W 


= o 1/3 r- 4/3 , 


v 2/3 -2/3 

X = a r , 


(A27) 

(A28) 

(A29) 

(A30) 


Y = a 1/3 r 2/3 , 

Z = -a- 1/3 r" 2/3 , 

or, explicitly, 

W = 3 1/s 2- 1/s (j + 1) 2/3 (1 -/3 PR )ao V“ 4/S x 

[2 1/2 / 3 i - 2~ 1,2 af 2 I(2 : 1)] “ 4/3 , (A31) 

X = 3 2/3 2- 2 / 3 (j + l) 4 / 3 (l-/3pR) 5 / 6 « 4/ V 2/3 x 

[2 1/2 / 3 i - 2~ 1,2 af 2 I(2 : 1)] “ 2/3 , (A32) 

y = 3 1 / 3 2- 1 / 3 (j + l) 2 / 3 (l- / 3 PR )- 1 / 6 Q 3/ V /3 x 

[2 1/2 / 3 i - 2“ 1/2 ao 2 I(2 : 1)] ^ , (A33) 

Z = -3- 1/3 2 1 / 8 (j^ + l)- 2/3 (l-^p R ) 1/6 a 0 - 3/ V 2/3 x 

1 — 2/3 


[2 1/2 / 3 i - 2“ 1/2 ao 2 I(2 : 1)]' 


(A34) 


where 1(2 : 1) = 1 if j = 1, and 0 otherwise. 

With J and b now defined, we can calculate the initial 
value of J, Jo, and the migration rate db/dt , which allows us 
to apply the capture probabilities, non-capture eccentricity 
kicks, and the distributions of initial libration widths, found 
by MW11. 

A further possible improvement over the MW11 model 
is the incorporation of eccentricity damping into the Hamil¬ 
tonian model. Under PR drag, e/e = 1.25a/a. Since the 
change in a when crossing a resonance is small, eccentricity 
does not decay greatly during the capture process; never¬ 
theless, we include this effect for completeness. Eccentricity 
damping is incorporated into the model by the addition of 
a term 


5 a dt 


2 a dt' ’ 


5 ( 

m p f 

21)9j V 


f M A 

- 2/3 

\m q J 



2/3 


= kbJ. 


(A35) 

(A36) 

(A37) 
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This introduces a third dimension to the parameter space. 
Rather than explore this thoroughly, we integrated two grids 
in Jo—db/dt space, one at k = 5 x 10 -5 and one at k — 0.02. 
These roughly correspond to the 2:1 resonance for m p \ = 
0.1 m© and 1000m© respectively. Despite the large difference 
in damping strengths, the capture probability is only weakly 
dependent on this (Fig |Al| . The eccentricity damping shifts 
the contours of the plot to the left, as some momentum 
J oc e 2 is lost prior to resonance passage. As we do not 
notice a significant effect, we do not add this complexity to 
our model. 


log IdB/dtl 


Resonant capture of dust 21 







